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ABSTRACT 

The Pluto system provides a unique local laboratory for the study of binaries with multiple low mass 
companions. In this paper, we study the orbital stability of P4, the most recently discovered moon 
in the Pluto system. This newfound companion orbits near the plane of the Pluto-Charon binary, 
roughly halfway between the two minor moons Nix and Hydra. We use a suite of few body integrations 
to constrain the masses of Nix and Hydra, and the orbital parameters of P4. For the system to remain 
stable over the age of the Solar System, the masses of Nix and Hydra likely do not exceed 5 x 10^^ kg 
and 9 x 10^^ kg, respectively. These upper limits assume a fixed mass ratio between Nix and Hydra at 
the value implied by their median optical brightness. Our study finds that stability is more sensitive 
to their total mass and that a downward revision of Charon's eccentricity (from our adopted value of 
0.0035) is unlikely to significantly affect our conclusions. Our upper limits are an order of magnitude 
below existing astrometric limits on the masses of Nix and Hydra. For a density at least that of ice, 
the albedos of Nix and Hydra would exceed 0.3. This constraint implies they are icy, as predicted by 
giant impact models. Even with these low masses, P4 only remains stable if its eccentricity e < 0.02. 
The 5:1 commensurability with Charon is particularly unstable. Combining stability constraints with 
the observed mean motion places the preferred orbit for P4 just exterior to the 5:1 resonance. These 
predictions will be tested when the New Horizons satellite visits Pluto. Based on the results for the 
Pluto-Charon system, we expect that circumbinary, multi-planet systems will be more widely spaced 
than their singleton counterparts. Further, circumbinary exoplanets close to the three-body stability 
boundary, such as those found by Kepler^ are less likely to have other companions nearby. 
Subject headings: planets and satellites: individual (Charon, Hydra, Nix, Pluto) — (stars:) planetary 
systems: formation — minor planets, asteroids — Kuiper belt — space vehicles 



1. INTRODUCTION 



New Horizons has one more objec t to study when 
it vis its Pluto in 2015 ( |Stern| |20Q8{ [Young fc Stern 
20To|. Hubble Space Telescope (HSTj imaging recently 



revealed a faint moon orbiting i n the p lane of the Pluto- 
Charon binary (S howalter et aT][2QTT] Sll). The moon. 



temporarily-named F4, orbits between the previously 
known o uter moons Nix and Hydra, which are 10 times 
brighter ([ Weaver et aL][2006[ ). 

Even before the discovery of P4, the Pluto system was 
dynamically intriguing. The orbital periods of Charon, 
Nix and Hydra are very close to a 1:4:6 ratio (Buie et al. 
[2006) . Despite thorough searching , no resonant lock has 
been identified ( [Tholenet al.|2008[ T08). The close prox- 
imity to resonance begs lor an explanation (as discussed 
further in the conclusions). The stakes are raised further 
by the proximity of P4 to a 5:1 commensurability with 
Charon (Sll). 

Efforts to understand these complex interactions are 
hindered by the difficulty of measuring the masses of Nix 
and Hydra. The tightest constraints come from the T08 
orbit solution, which fits astrometric data with four body 
integrations. However the low masses of Nix and Hydra 
make their dynamical perturbations difficult to measure. 
The T08 fits restrict Nix and Hydra to mass ratios below 
~ 10~^ with Pluto. These limits allow a wide range of 
plausible albedos and densities. 

The sizes of Nix and Hydra are indirectly constrained 
by their modest photometric variability, which suggests 



a roughly spherical shape (Stern et al. 2007[ S07). Nix 
and Hydra would have to be large, with diameters > 130 
km, for gravity to naturally make them spherical. Such 
large sizes are marginally inconsistent with the T08 mass 
constraints and quite inconsistent with our results. New 
Horizons should clarify lingering uncertainties about the 
size and shape of Nix and Hydra ^Young et al.[[2008 ). 

We exploit the precarious position of P4's orbit as an 
alternate approach to constraining system masses. The 
stability of P4 depends both on its uncertain orbital pa- 
rameters, and the strength of perturbations it receives 
from Nix and Hydra. We use a series of numerical in- 
tegrations to constrain simultaneously the masses of Nix 
and Hydra, and the orbit of P4. 
Prior to the announcement of P4, [Pires Dos Santos 



et al. (2011) investigated the stability of test particle or- 
bits in the Pluto-Charon system. Using the T08 masses, 
they found a pocket of low eccentricity orbits between 
Nix and Hydra that were stable for 10^ Charon periods 
(1800 years). Our paper assesses longer term stability 
(up to 3 X 10^ yr) for orbits in the vicinity of P4 for a 
range of Nix and Hydra masses. Ideally we would demon- 
strate stability up to the age of the Solar System. How- 
ever the short period of Charon (6.38 day) and the need 
to investigate many trial orbits make longer integrations 
costly. 

Because studying stability in the Pluto-Charon sys- 
tem places constraints on satellite masses, our dynami- 
cal studies also inform the composition of the bodies. By 
themselves, masses give reasonable constraints on surface 
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albedos and sizes. Once albedos and sizes are determined 
by New Horizons^ internal densities can be determined. 
The compositional diversity of bodies in the same system 



is in turn a powerful input to formation theories (|Benec- 
chi et al.|[2QQ9l |Stern 2009 ) 



As the prototypical — and also the second most mas- 
sive and best studied — Kuiper Belt object, Pluto rep- 
resents a critical stage in planet formation. As a transi- 
tional object, it hold clues both to the pro cesses that 
form the first pl anetesimals in gas disks (Chiang &" 
Youdin||2010t f^din 2010) and the collisional coagula- 



tion and destruction that in some pl aces produces planets 
and in others produces debris disks ([Kenyon fc Luu|1999 



Kenyon||2002[ [Kenyon fc Bromley|[2UTUi r 



We begm m ^ by describing parameters of the Pluto- 
Charon system, establishing which parameters are rea- 
sonably well constrained, and which must be varied in 
our simulations. Section |3| describes basic stability con- 
siderations for the Pluto-Charon system. Section [A[ finds 
the most circular orbits about the Pluto-Charon binary, a 
non-trivial task due to the non-Kepl erian nature of orbits 
about a binary (Lee & Peale 2006). Readers interested 
in the main results can skip to ^ which presents our 
numerical integrations of the Pluto system. In ^ we 
summarize our main findi ngs a nd discuss their implica- 
tions exoplanet studies in ^5.1 The appendix addresses 
technical aspects of the integrations including the ini- 
tialization of particles on the most circular orbit about a 
binary (^A| and details of the numerical code (^B|). 

2. PLUTO SYSTEM PARAMETERS 

2.1. Size and Mass of Nix, Hydra and P4 

Nix, Hydra and P4 are detected in reflected visible 
light. Their diameters depend on their unknown albedos, 
A, as 



^Nix 
^Hyd 



25A-^/^ km 
30A-^/^ km 
8^-^/2 km. 



(la) 
(lb) 
(Ic) 



To relate apparent magnitu des to size, we adopt a 
Charon-like phase coefficient ( Buie et aL[|1997[ T08). 

Mass estimates from photometry rely on an assumed 
density. With pi = p/{l g cm~^), the mass ratios rela- 
tive to Pluto are 



/iNix - 6.4 X 10-^1^"'/' 

/iHyd - 1.1 X 10-Vl^"'/' 

/iP4 ^ 2.0 X IQ-^piA- 



-3/2 



(2a) 
(2b) 
(2c) 



for Nix, Hydra, and P4. To express (without loss of 
generality) masses in terms of albedos, and not the more 

cumbersome Ap^ ' , we assume pi = 1, unless stated 
otherwise. Neglecting the possibility of high porosity, we 
also refer to A = pi = 1 as the "minimum mass" case. 

The orbit solution of T08 gives /iNix = (4.4±3.9) x 10"^ 
and /iHyd = (2.5 ± 4.9) x 10~^. The large uncertainties 
reflect the difficulty of measuring small mass ratios astro- 
metrically. Since the uncertainties extend towards or be- 
yond zerq^ the T08 solution places upper limits ~ 10~^ 

^ While the 1-cr errors on the mass of Nix do not quite extend 
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Fig. 1. — The dynamical environment of P4 is plotted in the 
space of semimajor axis, a, and eccentricity, e, relative to the Pluto- 
Charon center of mass. The green dot (with errorbars) gives the 
approximate location of P4. Vertical lines give nominal positions 
of mean motion resonances. Red and blue dashed lines give the 
nominal locations of the first order resonances with Hydra and Nix 
respectively. The black dot-dashed line denotes the 5:1 resonance 
with Charon (which nearly overlaps the 5:4 resonance with Nix). 
Black diagonal lines show the Nix and Hydra crossing trajectories, 
N-X and H-X, respectively. Magenta curves labelled Ct plot the 
critical Tisserand parameter relative to Nix and Hydra. The upper 
and lower Cx curves are for low and high masses of Nix and Hydra, 
i.e. high and low albedos as labelled. See text for details, 
on the mass ratios of Nix and Hydra. The equivalent 
albedo constraint is A > 0.04. Our stability calculations 
place much tighter constraints of A > 0.3 on Nix and 
Hydra. 

2.2. Orbit of P4 

The orbit of P4 is only loosely constrained. The discov- 
ery (Sll) announces P4 as "consistent with" a circular, 
coplanar orbit. Its mean motion from HST images span- 
ning ~ 20 days implies an orbital period of 32.1 ± 0.3 
days, corresponding to a period ratio of 5.03 ± 0.05 with 
Charon. 

For a Keplerian orbit about the Pluto- Charon barycen- 
ter (a good approximation at the 1% level), the measured 
period corresponds to a mean separation of 57400 ± 400 
km. Fig. [l] plots this orbit location as a dark green dot 
with error bars. The lighter green error bars show the 
larger range of orbits considered in our numerical simu- 
lations. 

The mean motion is consistent with the measured 
projected radial separation of 59000 ± 2000 km from 
Pluto (which is 2040 km from the barycenter). The pre- 
discovery HST images that date back to 2006 (Sll) as 
well as future observations from HST, ALMA and New 
Horizons will help constrain the orbit of P4. 

2.3. Orbital Ephemerides 



to zero mass, T08 discuss the difficulty of error estimation in their 
high dimensionality fits. 



Pluto's Circumbinary Chaos 



3 



TABLE 1 

Keplerian Orbital Elements^ and Cartesian State Vectors^ from Tholen et 

AL. (2008) 





Charon 


Nix 


Hydra 




Charon 


Nix 


Hydra 


Mass (kg)"" 


1.52e21 


(5.8el7) 


(3.2el7) 




0.1166 


(4.4e-5) 


(2.5e-5) 


a (km) 


19570.3 


49240 


65210 


X 


-0.1677 


-0.1744 


0.5202 


e 


0.0035 


0.0119 


0.0078 


y 


-0.2551 


0.7148 


-0.8822 


i 


96.168 


96.190 


96.362 


z 


0.0 


5.6384e-5 


3.4397e-3 




157.9 


244.3 


45.4 


X 


1.5995 


-1.0282 


1.0625 


M 


237.0 


129.80 


129.12 


y 


-1.0452 


-0.3580 


0.4458 


Q 


223.054 


223.202 


223.077 


z 


0.0 


-3.1683e-3 


-1.8596e-4 


P (days) 


6.3872 


25.49 


38.85 


p 
Ppc 


1 


3.99 


6.06 



^ Plutocentric for Charon and Barycentric for Nix and Hydra. The Nix and Hydra 
masses are the T08 values, not those used in this work. 
^ Cartesian state vectors are given in code units where G = Mp = Pj 
Pluto-Charon orbit has been rotated into the x-y plane. 



1. The 



Mp = 1.304e22 kg, /i is mass relative to Pluto. 



As summarized in Table [TJ we adopt the orbit solution 
of T08 for the orbits of Pluto, Charon, Nix and Hydra 
and the masses of Pluto and Charon. Ephemerides from 
JPL's HORIZONS systemp] are similar, but not identi- 
cal to the T08 solution. Orbit integrations started with 
the JPL (specifically the PLU021) ephemerides gave sta- 
tistically similar results to the T08 ephemerides. Thus, 
we only present results based the T08 orbit solution. 

Though the T08 (and also JPL) solutions assume spe- 
cific masses for Nix and Hydra, we vary the masses of 
Nix and Hydra without varying the orbit solution. The 
error introduced is hopefully modest compared to the as- 
trometric uncertainty. Future work could establish how 
orbit solutions vary with uncertain masses. 

Salient features of the system's orbital parameters, 
including circularity, coplanarity and the proximity of 
moons to resonances, are discussed in the introduction 
and evident in Table [T] Though low, the finite eccen- 
tricity of the Pluto-Charon orbit, ec = 3.5 x 10~^, has 
long been considered suspicious. Tidal interaction be- 
tween Plut o and Charon damp eg on short, ^ 10 Myr, 



timescales ( Dobrovolskis et al. 1997 Lithwick & Wu 
2008a). Nix and Hydra cannot excite ec above ~ 10~^ 
(LFOo); P4 can not contribute significantly either. Ex- 
ternal forcing — from solar and planetary tides, KBO 
fiybys and collision s — appear unable to explain the ob- 
served eccentricity ( Stern et al.||2003 ). 

After the submission of this manus cript, a lower ec - 
centricity for Charon was announced (Buie et al.||2012|). 
This new solution gives a 1-a limit of "3 km out of round" 
(Buie et al. 2 012), implying ec < 1.5 x lO""^ (at l-a). A 
preliminary investigation shows that setting ec = in 
the T08 orbit solution has a minor effect on orbital sta- 



bility (^4.4). 



3. BASIC STABILITY CONSIDERATIONS 

This section summarizes previous results on orbital 
stability that are relevant to the Pluto system in par- 
ticular and to multi-satellite systems about a binary in 
general. While these general considerations cannot pre- 
dict the stability of P4, they are useful in understanding 
the dynamical environment of P4 as shown in Fig. [l] 
and in interpreting the numerical simulations in Q Our 
most significant general conclusion concerns the compar- 

^ Accessible via http://ssd.jpl.nasa.gov/horizons.cgi 



ison of previous work on multi-planet stability about a 
single star (^3.3) to our simulation of the Pluto-Charon 
circumbinary system. The comparison shows that binary 
perturbations can significantly enhance the destabilizing 
effects of satellite interactions, which is a small step to a 
more complete understanding of these complex dynami- 
cal interactions. 

3.1. Stability of Circumbinary Orbits 
The stability of a test particle orbiting an inner bi- 



nary is well-studied ( [Szebehely 



1980 Dvorak 1986), in- 



cluding the system atic numerical investigation of Hol man| 
fc WiegertI ( |1999[ hereafter HW). The HW results show 
that INix, Hydra and P4 are all individually stable about 
the Pluto-Charon binary. Based on the results of HW 
(see their equation 3 and/or their table 7), period ra- 
tios below 2.8 are unstable in the Pluto-Charon system. 
Nix's period ratio of ~ 4 is 41% larger than the stabil- 
ity boundary, and also avoids a possible instability strip 
near the 3 : 1 resonance. 

3.2. Stability of Two Satellites 

We now temporarily ignore the important fact that 
Pluto and Charon are a binary and consider the 3- 
body stability of two satellites about a central mass that 
places the combined mass of Pluto and Charon at their 
bary center. This approximate problem allows us to sim- 
ply demonstrate how strong interactions with Nix or Hy- 
dra restrict the allowed orbits of P4, as in Fig.[l] and also 
aids the interpretation of numerical results. 

The well known stability criterion for initially circular 
orbits with semimajor axes ai and a2 = ai + A, re- 
quires A > 3.5i?H in the restricted thre e body probleni 



when one of the satellites is massless (iGladman 1993 



which also analyzes two massive satellites). That is, two 
satellites on circular orbits are stable if their orbits are 
separated by more than 3.5 Hill radii. 



(3) 



of the massive satellite. 

The separation of P4 from Nix or Hydra in terms of Hill 
radii depends on the assumed masses for Nix and Hydra. 
For the mass scalings in Equation ([2| , P4 is separated 
from Nix by 29^/Am^ and from Hydra by 17^/A^^ Hill 
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radii (of Nix and Hydra respectively). Separation by 
fewer than S.bRn would require either Anix < 0.014 or 
^Hyd < 0.042. Such low albedos and thus high masses 
are inconsistent with the T08 (la) upper mass limits. 
The more stringent mass limits that we find ensure that 
P4 is stable by the Hill stability criterion for circular 
orbits. 

The stability of an eccentric test mass follows from the 
(approximately) conserved Tisserand parameter 



a 

2a 



(4) 



where a, e, i (inclination) refer to the test body (P4) and 
a' refers to the perturber (Nix or Hydra). The stability 
boundary of circular, coplanar (e = z = 0) orbits at 
\a — a'\ = 3.5i?H defines a critical Ct- (Restricted 3 
body) stability for arbitrary e and i holds for Ct above 
the critical value. 

Fig. [l] plots the critical Ct curves, sometimes called 
"Tisserand tails." The outer and inner tails of Nix and 
Hydra, respectively, are relevant for interactions with P4. 
For low Nix and Hydra masses (the A = 1 curves) the 
tails would only cross an eccentric P4, ep4 ^0.1. With 
higher Nix and Hydra masses (the A = 0.05 curves). 
Hydra's Tisserand tail intersects plausible P4 orbits at 
lower ec, especially if P4 lies on the outer edge of allowed 
orbits. While a considerable simplification, these consid- 
erations from the restricted three body problem demon- 
strate how higher masses for Nix and Hydra severely re- 
strict the allowed orbits of P4. 

3.3. Stability of Three Satellites 

Still ignoring perturbations from the central binary, we 
now consider three interacting satellites about a central 
mass. Sharp stability boundaries no longer exist, but the 
timescale to orbit crossing can be studied numerically. 

Most investigations of multi-planet stability consider 
roughly equal mass compa nions with evenly spaced or- 



bits, in terms of Hill radii. [Chambers et al.| (p^96j here- 



after C96) defined the mutual Hill radius, as 



1/3 



CLi + CL, 



+ 1 



(5) 



for neighboring planets {i and i + 1). This definition does 
not reduce to the standard Hill radius when /i^ ^ 0. 
This deficiency is readily corrected by a mass-weighted 
averaging of the semi-major axes. This distinction is 
insignificant when the satellite masses are similar, but is 
relevant here due to the low mass of P4. 

C96 measured the orbit crossing timescale, tc, for sys- 
tems of three equal mass planets, with planet-star mass 
ratios ranging from /i = 10~^ to /i = 10~^. Relative 
to the period of the inner planet. Pi, we fit the data in 
Figure 4 of C96 as 

logio(tc/Pi) = -9.11 + 4.39AV1/12 _ i.071og(/i) , (6) 

where A' is the orbit separation in mut ual Hill radii and 
the fu nctional form is motivated by the |Faber fc Quillen 
(2007) study of systems of > 3 planets. The mass scaling 



A'/i-^/-^^ cx /i differs from t he /i Hill ra dius scaling 



Unequal masses (of P4 in particular) make direct ap- 
plication of these results difficult. To allow simple esti- 
mates, we make a range of possible assumptions: /i as the 
average of all 3 satellites or of just Nix and Hydra; the 
mutual Hill radius defined as in Equation ([5| or with den- 
sity weighted a^. In all cases the minimum mass A = 1 
case is stable, with crossing times > 10^^ yrs. For the 
higher mass A = 0.05 case, tc ^ 3 x 10^ yrs, implying 
rapid instability. 

Even neglecting binary perturbations, the stability of 
Pluto's minor moons depends strongly on their assumed 
masses. Binary perturbations accelerate the destabiliza- 
tion. We show in section |4] that for the high mass case 
with A = 0.05, the simulated crossing time for P4 is 
< 10^ yrs, over a thousand times faster than the above 
estimate neglecting the central binary. 

3.4. External Perturbations & Tides 

In this paper we approximate the Pluto system as a 
set of isolated point masses. Since Pluto's Hill sphere ex- 
tends '^120 times further than Hydra's orbit, solar tides 
are a minor effect. The protective 3:2 resonance between 
Neptune and Pluto helps to weaken the dominant plan- 
etary perturbation. Nevertheless, followup work should 
test the effect of weak external perturbations on long 
term stability. 

Collisions with interloping KBOs could significantly 
perturb the weakly bound outer moons. If the Kuiper 
belt was massive in its youth, c ollisions would unbind 



100 km class binary KBOs (Nesvorny et al. 2011 



Kavelaars[2012 ). Pluto's moons, F4 in particu- 



Parker 

lar, could be destabilized by collisions that only modestly 
perturb its eccentricity. Combining the dynamical stabil- 
ity of Pluto's moons with collisional perturbations could 
be a powerful constraint on both the Pluto system and 
the collisional environment of the Kuiper belt. 

Tides in the Pluto system are dominated by interac- 
tions between Pluto and Charon, which are locked in a 
dual synchronous s tate, with both spin perio ds equalling 



the orbital period ( [Dobrovolskis et al. p!997 ). In princi- 
ple, the eccentricities of minor moons snould be damped 
by exciting the eccentricity of Charon, which then suf- 
fers tidal dissipation. However, iLithwick & Wu (2008b) 



due to three-body resonances ( |Quillen||2QTT] ). 



conclude that this damping mechanism is weak. 

4. RESULTS FOR LONG-TERM STABILITY 
4.1. 4+^ Body Integrations 

We study the stability of P4 with a suite of 4 + N body 
integrations. In these simulations, the four massive bod- 
ies are Pluto, Charon, Nix and Hydra. Treating P4 as 
a massless test particle allows simultaneous investigation 
of many trial orbits. We follow each test particle until it 
crosses the orbit of either Nix or Hydra. (No significant 
perturbations to the orbits of the massive bodies occurs 
in the simulations.) Integrations were performed using 
the 15th order Radau integrator (Everhart 1985), as im- 
plemented by the Swifter software packager] Details con- 
cerning code performance are deferred to the appendix. 

Between different sets of simulations, we vary the un- 
certain masses of Nix and Hydra, keeping their mass ratio 
fixed at MNix/^Hyd = 0.575, the value implied by their 

^ Publicly available at http://www.boulder.swri.edu/swifter/. 
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Fig. 2. — Median crossing time, tc, versus initial a [as Aa : 



(5 X 10^ km)] and e of P4. Tlie mass of Nix and Hydra is labelled as 



both A (their shared albedo for p = 1 g cm~^) and mi (mass relative to the minumum mass at A = 1). The color scale for the tc is above 
each panel. Crossing times are longer for lower masses of Nix and Hydra and lower eccentricities of P4. The semi-major axis dependence 
is more complicated, due to the effect of resonances. See text for details. 



median optical brightness. This mass ratio assumes that 
Nix and Hydra share the same density and albedo, which 
should be expected for standard formation scenarios]^ 
The range of Nix and Hydra masses considered corre- 
sponds (for a density of 1 g cm~^) to albedos from 0.03 
to 0.24. Integrations were run with albedos up to 0.4 
(/iNix ^ 10"^), however only 10 - 30% of test particles 
suffered orbit crossing after ~ 10^ Charon orbits, making 
systematic investigations too expensive. 

The initial conditions for the massive bodies are given 
in Table [l] The initial orbits for the test particles were 
chosen by two differen t me thods. The first suite of sim- 
ulations, described in §4.2[ populated P4 orbits by ran- 
domly sampling Keplerian ele ment s. The second suite 
of simulations, summarized in §4.3| initializes P4 on the 
"most-circular" orbits about the Pluto- Charon binary. 



4.2. Uniformly Sampled Keplerian Orbits 

For each adopted mass of Nix and Hydra, we integrate 
5000 P4 orbits with initial conditions chosen randomly 
from a uniform distribution of Keplerian osculating ele- 



While the best fit T08 masses have a different ratio, large 
uncertainties mean the ratio has little significance. Furthermore, 
integrations showed that the T08 masses destabilize P4 orbits on 
a median timescale of < 10^ years. 



mentsj^ The semimajor axes are restricted to the range 

56632 km < a < 58832 km , (7) 

or a/apc = 2.89 - 3.01. Eccentricity and inclination 
restricted to e < 0.05 and i < 0.5°. The modest range in i 
was insignificant for our results and will not be discussed 
further. Other Keplerian angles (argument of pericenter, 
longitude of ascending node and mean longitude) were 
sampled over the full (27r) range. Fig. [l] compares the 
sampled orbits to the nominal location or mean motion 
resonances as well as Nix and Hydra crossing trajectories. 

4.2.1. Mapping a — e Space 

Fig. [2] maps the median stability timescale versus ini- 
tial a and e for several Nix and Hydra masses. As shown 
by the colorbars on each map, crossing times increase 
significantly as the mass of Nix and Hydra drops from 
high to low (upper left to lower right plots). Crossing 
times are always short at higher a and e (upper right of 
each plot). Interactions with Hydra are the likely culprit, 
consistent with the Tisserand parameter curves in Fig. [l] 

At low eccentricity, crossing times are longest, and dis- 
play a complex dependence on a. Resonances play a key 
role. Indee d resonance overlap is a generic cause of or- 
bital chaos ( Mudryk fc Wu|2006 ) . From the approximate 



^ We use Jacobian elements are measured relative to the barycen- 
ter of Pluto, Charon and Nix. 
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Fig. 3. — Crossing timescale for P4 versus the masses of Nix 
and Hydra. The bottom axis scales masses relative to the albedo 
one case. The top axis shows the corresponding albedo (for 
p = 1 g cm~^). Circles give the crossing times from simulations. 
Dashed lines are power law fits to the mass dependence. Colors 
denote different sets of initial P4 orbits. From bottom to top, me- 
dian crossing times are shown for: (blue:) the full range of orbital 
parameters; (green:) e < 0.015; (black) a stable "core" with both 
e < 0.015 and 5.70 x 10^ km < a < 5.75 x 10^ km. Finally red data 
points give the 90th percentile of longest lived orbits in the core. 
If extrapolation can be trusted, Nix and Hydra require A > 0.5 to 
ensure the stability of P4 over the age of the Solar System. 

resonance locations in Fig. [l] specific resonances can be 
implicated. 

The 5:1 resonance with Charon helps destabilize the 
region near a ^ 5.8 x 10^ km (Aa ~ 8.0[xl0^ km] 
as labelled in Fig. [2| at lower masses. The 4:5 reso- 
nance with Hydra shortens crossing times between 5.65 < 
a/{10^ km) < 5.70 (i.e. 6.5 < Aa/(10^ km) < 7.0), for 
the highest mass (upper left). This change in the promi- 
nence of different resonances is expected as the masses 
of Nix and Hydra vary relative to Pluto- Charon. 

4.2.2. Median Crossing Times: Measured and Extrapolated 

Fig. [3] shows how crossing times scale with the masses 
of Nix and Hydra. Median timescales are plotted for 
three subsets of the initial orbital parameters: (1) all 
initial orbits, (2) only e < 0.015, and (3) the stable "core" 
of parameters with both e < 0.015 and 5.70 x 10^ km < 
a < 5.75 X 10^ km. For the "core" sample, we also plot 
the 90th percentile of crossing time (beyond which only 
10% of particles survive). At any mass, crossing times 
increase as cuts become more selective. The "missing" 
data points for low mass cases occur where P4 orbits were 
so stable that the median (or 90th percentile) timescale 
was not reached after 10^ Pluto-Charon periods. 

The dependence of crossing time on the mass, m, of 
Nix and Hydra is well described by a powerlaw 

tc (X . (8) 

The best-fit indices are 7 ^ —3.6, —4.3, —4.4, —4.6 for 
the bottom to top curves in Fig. |3] The larger scatter 



about the 90th percentile "core" powerlaw is partly due 
to Poisson noise, as fewer orbits contribute to this re- 
strictive sample. Physical scatter caused by the shifting 
locations of the most stable regions could also contribute. 

No known theoretical explanation exists for such a 
powerlaw scaling. Indeed, if crossing times scale with 
Hill radius as in Equation (pi) , the mass dependence is not 
a simple powerlaw — instead the local powerlaw index 



would steepen towards lower masses. However, [Duncan 
&: Lissauer[ ( |199 7) found that orbit crossing timescales for 



the Uranian moons also exhibit a powerlaw dependence 
on satellite mass. 

Extrapolation along these powerlaws allows us to spec- 
ulate about the stability of P4 for very low masses of 
Nix and Hydra. While such extrapolation is inherently 
uncertain, low masses are much more costly to simu- 
late directly due to longer crossing timescales. The goal 
of extrapolation is to estimate the "allowed" masses of 
Nix and Hydra, where allowed means that P4 crossing 
timescales exceed the age of the Solar System. 

Including all sampled P4 orbits (the blue curve in 
Fig. [3]), the allowed masses of Nix and Hydra are implau- 
siblylow — below the minimum set by A = 1. However 
for low eccentricity P4 orbits (all other curves), a range 
of plausible Nix and Hydra masses are allowed. From 
extrapolation of the 90th percentile "core" sample (red 
curve), the lowest allowed albedo is A ~ 0.5. Lower albe- 
dos are possible for special orbits, including the "most 
circular" orbits that we present next. 

4.3. Stability of the ''Most Circular'' Orbits 

The greater stability of low eccentricity P4 orbits 
shown above motivates the use of the "most circular" 
orbits about Pluto-Charon as an initial condition. These 
special orbits are not simply found by setting the Kep- 
lerian e = 0. Appendix [A] describes techniques to find 
these orbits. 

We integrated trial orbits for 48 distinct P4 periods: 
41 are uniformly spaced with period ratios (relative to 
Charon) from 4.79 to 5.31 and the remaining 7 provide 
more dense sampling near the 5:1 resonance. At each 
orbital period 25 trial orbits where integrated. These 25 
orbits are shifted in orbital phase along the most circular 



orbit by one Charon period each, as described in ^ A.2 

Fig. [4| shows the crossing times of the most circular 
orbits, for all orbital periods and for different Nix and 
Hydra masses. The median crossing times at each period 
are given by square symbols, with errorbars spanning 
the 10th and 90th percentile. Lower limits are given 
where simulations did not run long enough to determine 
a timescale. 

Crossing times become longer as mass decreases, and 
the dependance on orbital period (or a) is quite complex. 
This general behavior agrees with the Keplerian initial 
conditions ( ^4.2| ). Fig. [1] shows further that the period 
dependence is more varied and irregular for lower masses 
of Nix and Hydra. While it is difficult to understand all 
these fluctuations in detail, the trend is intuitively con- 
sistent with weaker resonant forcing being more narrowly 
focused at specific orbital periods. 

Fig. |4] shows a dramatic drop in crossing times near 
Charon's 5:1 resonance (C5:l). As also seen in Fig. [2 
this instability strip appears for lower masses of Nix anc. 
Hydra. At higher masses. Nix and Hydra perturbations 
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Fig. 4. — Crossing times versus orbital period for P4 when ini- 
tialized on the "most circular" orbits about Pluto-Charon. The 
masses of Nix and Hydra (taking the same values as Figs. [2] and [sj 
decrease from the bottom to top sets of colored symbols. Ivlarkers 
with errorbars denote the median timescale and the 10th - 90th 
percentile. Triangles indicate lower limits. Nominal locations of 
mean motion resonances with Charon, Nix and Hydra are shown 
along the bottom axis. For lower masses, a narrow instability strip 
exists at the 5:1 commensurability with Charon. The grey hori- 
zontal errorbar shows the observed mean motion of P4 (Sll). 



wash out the influence of C5:l. Longer- hved P4 orbits ex- 
ist both outside and inside the narrow instabihty strip at 
C5:l Within the observations constraints, shown by the 
horizontal error bar, our analysis favors locations outside 
C5:l. 

Fig. [5] plots crossing times versus the mass of Nix and 
Hydra, with powerlaw fits overplotted. We restrict the 
range of periods to the observational constraint (Sll), 
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Fig. 5. — Similar to Fig. [3] but for crossing times, tc, of the "most 
circular" initial orbits witn period ratios (to Charon) between 4.93 
and 5.13. Blue squares give the median over both period and 
initial phase. Green squares (and red squares) give the longest 
of the median (and 90th percentile) tc at each period. Yellow 
squares give the longest tc of all orbits. For reference the solid 
line gives powerlaw fit to the median tc for the "core" sample of 
Keplerian initial conditions. The most circular orbits give longer 
tc and thereby allow larger Nix and Hydra masses. 



but with double the uncertainty for inclusiveness. The 
statistical measures of tc include: (1) "median," which 
takes the median of the values given by the symbols in 
Fig. [4] (themselves median values over phases); (2) "max 
median," the longest phase- median tc, i-e. highest sym- 
bol; (3) "max 90%," the longest 90th percentile tc, i.e. 
highest upper errorbar; and finally (4) "longest," simply 
the longest tc at any phase or period considered. 

For reference, the median tc for the relatively stable 
"core" sample of initial Keplerian parameters is over- 
plotted. Compared to this reference case, the cross- 
ing timescales for the most circular orbits are signifi- 
cantly longer, especially toward the (more realistic) lower 
masses of Nix and Hydra. Longer crossing times equate 
to higher allowed masses (and lower albedos) for Nix and 
Hydra. 

Extrapolating along the powerlaw fits in Fig. 5] shows 
that A > 0.3 is needed to achieve tc > 4 Gyr. This limit 
is more inclusive than the A ^ 0.5 found in ^4.2.2[ We 
cannot definitely rule out even lower albedos. As already 
discussed, extrapolation could prove misleading. 

Characterizing the most stable orbit is difficult. We 
do not base our estimate on the absolute longest lived 
orbits, which loosen our constraints. As shown in Fig. [5j 
the longest tc's do not follow a simple powerlaw and are 
therefore unreliable for extrapolation. Even without ex- 
trapolation, the longest tc's are highly subject to sam- 
pling, especially considering the pronounced period and 
phase dependence shown in Fig. [4j It is also unclear if 
P4 is likely to inhabit the most stable orbits, especially if 
those orbits occupy a tiny volume of pha se sp ace. Small 
neglected effects (such as collisions, see ^3.4) could eas- 
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ily remove P4 from narrow pockets of parameter space. 
Thus, we base our constraints not on the absolutely most 
stable orbit, but on an average of orbits among the most 
stable. 

4.4. Parameter Tests: Satellite Mass Ratios and Binary 
Eccentricity 

We include a preliminary investigation of the effect of 
Nix and Hydra's mass ratio on the orbital stability of P4. 
Our main set of simulations fixed MNix/^Hyd = 0.575. 
Here we also investigate MNix/^Hyd = 0.3 and 1.2, i.e. 
roughly half and double the fiducial value. For each mass 
ratio, we consider 4 values of the fixed total mass equal to 
the 4 highest mass cases in our main set of simulations. 
As in ^4.3 , we initialize the P4 orbits on the most circular 
orbits, with the same sampling of orbital periods and 
phases. 

Fig. |6] plots P4 lifetimes for these runs with different 
mass ratios. The dependence of orbital stability on the 
Nix/Hydra mass ratio is relatively minor compared to 
the stronger dependancies on Nix and Hydra's total mass 
and P4's orbital period. The decrease in orbital stabil- 
ity near the 5:1 commensurability - probably the most 
relevant feature - is evident in the lower mass runs for 
all mass ratios investigated, as are several other trends 
with period. The most prominent effect of the mass ra- 
tio is an increase in P4 crossing times at longer periods 
( > 5.1Pcharon) for larger mass ratios, i.e. smaller Hydra 
masses. This behavior is unsurprising, and is qualita- 
tively explained by 3 body interactions. As shown in 
Fig. [l] a more massive Hydra strongly perturbs orbits 
exterior to P4's measured orbit. 

For the two lowest total mass runs, the variation in 
median crossing timescales was < 50% between the dif- 
ferent mass ratios, compared to the fiducial case. Our 
preliminary investigation suggests that it is difficult for 
stability studies to precisely constrain Nix and Hydra's 
mass ratio at present. However, ever-improving orbit 
determinations should help, both on their own and in 
combination with stability studies. 

We also investigated the effect of setting Charon's ec- 
centrici ty to zero. We conducted this test with the meth- 



ods of ^.2, i.e. uniform sampling of Keplerian osculat- 
ing elements. As above, we investigated the four highest 
mass satellite cases. We find that reducing ec to zero has 
a modest stabilizing effect. The median crossing time in- 
creased by < 50% (40%, 30%, 17% and 42% for the high- 
est to lowest satellite masses considered). Future work 
should aim to develop a more systematic understanding 
of how binary eccentricity, and other orbital parameters, 
affect stability. 

5. SUMMARY AND DISCUSSION 

We study the long term stability of P4, the temporary 
name for the moon orbiting the Pluto-Charon binary be- 
tween Nix and Hydra (Sll). Our numerical integrations 
constrain both the orbit of P4 and the masses of Nix and 
Hydra. These constraints are coupled, so improved de- 
termination of P4's orbit will help refine the masses of 
Nix and Hydra and vice-versa. We summarize our main 
results: 

• Low eccentricity orbits of P4 are significantly more 
stable. Our integrations strongly disfavor P4 orbits 
with e > 0.02, as shown in Figs. [2] and |3] 



• Period ratios (of P4 to Charon) between 4.98 and 
5.01 are unstable on short timescales. Slightly 
larger or smaller period ratios are significantly 
more stable, as shown in Fig. [4j Combined with 
the observed mean motion (Sllj, our results favor 
orbits just outside the 5:1 commensurability with 
Charon. 

• Even the most stable P4 orbits only survive if Nix 
and Hydra are sufficiently low in mass. We esti- 
mate MNix < 5 X 10^^ kg and Muyd < 9 x 10^^ kg 
are required for the stability of P4 over the age of 
the Solar System. This constraint holds the mass 
ratio between Nix and Hydra fixed at the value im- 
plied by their mean brightness. 

• The albedos of Nix and Hydra are correspondingly 
constrained to A > 0.3, assuming an internal den- 
sity of 1 g cm~^. Higher density rocky bodies would 
require even higher albedos. 

• The above mass and albedo constraints rely on ex- 
trapolation of simulations with higher masses, as 
shown in Fig. [5] Direct simulations alone disfavor 
A < 0.16 for which the orbit crossing time of P4 is 

< 10' yr. 

Our mass limits based on the stability of P4 are a fac- 
tor of 20 and 10 lower than the (1-cr) astrometric upper 
limits of T08. The rendezvous of the New Horizons satel- 
lite with the Pluto syst em in July 2015 a s well as ongo- 



ing Hubble astrometry (Buie et al 

astrometric mass constraints' 



improve 



2012|) should greatly 
Neglecting P4, 

Beauvalet et al.| (l2012| combine current data with simu- 
lated New Horizons observations to show that mass er- 
rors on Hydra will be reduced to ~ 4 x 10^^ kg. This limit 
is already small enough to test our predictions. Hopefully 
the inclusion of P4 will further tighten astrometric mass 
constraints. Ultimately, combining astrometry with long 
term stability should provide the tightest and most ro- 
bust dynamical constraints. 

Our results generally support the leading model for the 
origin of the Pluto system : a fflant inipact that produces 
the Pluto-Charon binary ( |McKinnon||1989| |Canup 20 Q5[) 
and t he debris that forms its coplanar moons (Stern e t al.| 
2006). The low eccentricity of P4 requires collisional 
damping. The inferred high albedo of Nix and Hydra 
is consistent with these being icy bodies. In the model 
of jCanup (2011), the collision of differentiated bodies 
(or rock with icy mantles) forms Pluto and Charon plus 
a pure ice debris disk from which the moons can accu- 
mulate. Collisional stripping of an icy mantle similarly 
explains many properties of the dwarf plan et Haumea: 
rapid rotation, collisional family members (Bro wn et al. 
2007[), and two high albedo icy moons ( |Ragozzine 
Brown||2U09|) . A possible weakness of the collisional see- 



nario is that the debris disk forms much closer to PlutO' 
Charon than Nix's current orbit. Outward orbital migra- 
tion has been proposed, but issues regarding simultane- 
ous migration of multiple moons — whic h are likely more 
severe with P4 — remain unresolved (Ward fc Canup 



2006[ [Lithwick fc Wu|[2Q08a[ |Canup||20lT| r 



The capture scenario for the system is implausible, es- 
pecially for the circular and coplanar minor moons. How- 
ever, the gravitational collapse scenario for the formation 
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to-pr imary mass ratios of 0.3, 0.97 and 0.91 (Welsh et al. 
2012). The period ratios of the planets to the inner bi- 
naries are 5.6, 10.4 and 6.3 (respectively). Unlike the 
Pluto system, there is no preference for being near n:l 
resonances. However the Kepler multiples around single 
stars do show a preference (as yet unexpl ained) for bein^ 



near low order mean motion resonances (|Fabrycky et al. 
20121). 



Fig. 6. — Similar to Fig.JJl but varying the mass ratio of Nix to Hydra (0.3,0.575, 1.2 for left, center and right panels) for four values 
of the combined Nix and Hydra masses (with the same values and colorscale as Fig. pi). The qualitative stability behavior is similar, 
including the decreased stability near 5:1 with Charon for the two lowest total masses (cyan and yellow data). At large periods, stability 
times increase for higher mass ratios, i.e. relatively smaller Hydra masses. 

of lower mass Kuiper belt binaries mi^ht also ap ply to 
the Pluto-Charon system (Nesvorny et al. 2010). This 
binary formation model is a natural extension or a lead- 
ing planetesimal formation theory: solid debris accumu- 
lates via the streaming instability (Youdin & Goodman 
2005 Youdin fc Johansen||20^7|) and subs equently un- 
der^Qes~jsi:ravitation al collapse ( [Youdin fc Shu^,2002] Jo- 
' hansen et al.||2009 ). Rotational fission is the extra twist 
needed tor binaries. A particularly massive clump would 
be required for the formation of Pluto and Charon. How- 
ever massive clumps form by merging in streaming insta- 
bility simulations (Johansen fc Youdin^2007) and massive 
particle ring s could form via related secular gravit ational 
instabilities ( |Youdin||2011a| [Shariff fc Cuzzifc oiT). 

To explain F lut o ' s outer moons, some of tne collapsing 
material must accumulate in a disk around the binary. 
Indeed a possible benefit of the scenario is that higher 
angular momentum collapsing material could more read- 
ily produce distant moons. While plausible, this specific 
scenario has not been modeled in detail. 



TTep/er's circumbinary planets are relatively close to 
the circumbinary stability boundary (HW). The period 
ratio of Keper 16b is only 14 % beyong this limit fW elshl 
et al.|2012[ ) compared to 41% for Nix. Of course, the ob- 
servational bias es of transit sur veys strongly favor shorter 
period planets (|Youdin||2011b|). 



5.1. Pluto- Charon as Exoplanet Host Stars 

As described in ^ the dynamics of the Pluto-Charon 
system are broadly relevant to our understanding of cir- 
cumbinary dynamics. We conclude with a brief compar- 
ison to the circumbinary planets discovered by Kepler^ 
and emphasize that Pluto is a guide to the circumbinary 
multi-planet systems that have yet to be discovered. 

To make an analogy with planetary systems, we may 
scale the mass of Pluto to a Solar mass. Charon then 
equates to a 0.1 2M0 red dwarf. Nix, Hydra and P4 
would then have scaled minimum masses of ^ 2 MMars, 
~ 3.5 MMars and ~ 3 Mp respectively. Kepler 16, 34 
and 35 contain planets with masses of 1.1, 0.7 and 0.4 
times Saturn (respectively) around binaries of secondary- 



Kepler has revealed that compac t multi-planet systein s 
are common around single stars (Lissauer et al. |2011 ). 
However for circumbinary systems, a planet too close to 
the stability boundary is unlikely to have other planets 
nearby. Thus circumbinary planets detected at a safer 
distance from the stability boundary are more likely to be 
in circumbinary multi-planet systems — the true Pluto 
analogs. 

Our study shows that multi-planet scattering is more 
violent around a binary than a single star. Thus orbital 
stability will require larger spacing between circumbinary 
planets, making detection of multi-planet systems more 
challenging. Nevertheless, Pluto's rich set of companions 
bodes well for a fruitful search. 
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APPENDIX 



A. "MOST CIRCULAR" ORBITS ABOUT A BINARY 
Section [43] presents integrations of the Pluto system with P4 initially placed on the most-circular orbit about the 



Pluto-Charon binary. Lee & Peale (2006) develop an analytic theory for the orbits of test masses about Pluto and 



Charon (or more generally any circular binary). Their theory separates orbital motion into three components: (1) 
circular motion of the guiding center about the barycenter, (2) oscillations forced by binary motion, and (3) the free 
or epicyclic eccentricity. Setting the epicyclic eccentricity (distinct from Keplerian eccentricity) to zero gives the most 
circular orbits. Their theory applies to a circular binary orbit, for which the potential is constant in a rotating frame. 
To account for the small (but probably overes timated) ecce ntricity of Charon in the T08 orbit solution, we instead use 



the invariant loop method of Pichardo et al.| (2005), which Lithwick & Wu (|2008a ) applied to the Pluto system. We 
demonstrate below that the adopted value of Charon's eccentricity gives only a minor correction to the |Lee fc Peale 
( [2006, ) solutions for the most circular orbits. 



A.l. Finding Invariant Loops 

In time-periodic potentials — such as an eccentric binary — special orbits form a closed loop when plotted strobo- 
scopically, i.e. when a snapshot is taken once per orbit of the potential. For concreteness we take this snapshot when 
the binary orbit is at periapse. These special, stroboscopically closed orbits are the most circular orbits. In the case 
of resonant orbits, the stroboscopically closed loop breaks into a chain of islands. 

Invariant loops in the plane of the binary can be found by iterative methods. Integrations start when the binary is 
at periapse and the test particle is aligned with the binary. By symmetry, the radial velocity is zero at this location for 
an invariant loop. Thus at a given radial distance, the only initial condition to vary is the azimuthal velocity, Vq- For 
the correct choice of ^o, the stroboscopic orbit forms a ID curve. To iterate towards this orbit we vary Vo to minimize 
the radi al dispersi on of the stroboscopic orbit. 



[Lithwick fc Wu (|2008a) measured the radial dispersion in narrow azimuthal bins so that the invariant loop has a 
nearly constant radial distance in the bin. We found faster and more reliable convergence by fitting the stroboscopic 
orbits to a (fourth order) cosine series, and measuring the radial dispersion about that best fit orbit. 



A. 2. Implementation as Initial Conditions 

In §4.3[ we present integrations where Nix and Hydra perturb P4 from these most circular orbits. The initial orbital 
phase along the most circular orbit must be consistent with the initial orbital phase of the binary. This initial phase 
is not unique as P4 can be advanced along the orbit by an integer number of Charon periods. For different orbital 
phases, the position of P4 changes relative to Nix and Hydra. Due to the sensitive dependence on initial conditions, 
stability timescales also change. Thus our simulations sample many (typically 25) initial phases of P4 for every trial 
orbital period. 



A. 3. Properties of the Most Circular Orbits 

Fig. [7| plots the time evolution of the most circular test particle orbits in the vicinity of P4. These orbits only include 
the perturbations from the Pluto- Charon binary, and not those due to Nix or Hydra. Each curve represents a different 
orbital period, ranging from 4.8 to 5.2 Pluto-Charon periods. None of these orbits are in a 5:1 resonan ce with Charon, 
and the effect of the 5:1 is negligible due to the low ec, unlike the study of the 4:1 for a high ec by |Lithwick fc Wu[ 
( [2QQ8ap . 

'I'he left panel plots radial distance from the Pluto-Charon barycenter. The non-circularity, given by the magnitude 
of the radial excursions relative to the mean r, is Ar/(2r) ~ 2 x 10~^ . The periodic radial oscillations have well- 
understood timescales (Lee & Peale 2006). The faster oscillations correspond to the synodic period with Charon 
(~ 5/4 Charon periods) and the first harmonic at half that period. The slower oscillations are on the epicyclic period 
(~ 5 Charon periods) of the test p article itself. The m ost circular orbits about a circular binary do not have these 
longer period epicyclic oscillations (Lee & Peale 2006). The epicyclic oscillations are only modestly larger than the 
synodic timescale variations. When Nix and Hydra are introduced to numerical integrations they rapidly excite a 
larger epicyclic eccentricity. Thus Charon's adopted eccentricity is not expected to strongly affect orbital stability. 
The integrations in ^4.4 affirm this expectation. 

The middle and right panels of Fig. [t] plot the osculating (instantaneous) Keplerian a and e, respectively, about the 
barycenter to demonstrate the non-KepIerian nature of these orbits. The Keplerian a is both systematically larger than 
the cylindrical radius and varies more significantly with time. The Keplerian e oscillates and reaches values ~ 0.015, 
much larger than the actual radial excursions. This plot demonstrates clearly that setting the osculating e = does 
not give the most circular orbit, and why e had to be increased above ~ 0.01 in Fig. [2] to noticeably affect orbital 
stability. 
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Fig. 7. — Evolution of the "most circular" test particle orbits about Pluto-Charon versus time. The orbits are in the vicinity of P4 and 
neglect perturbations from Nix and Hydra. Left: Radial distance from the Pluto-Charon barycenter. The denser packing of lines near the 
middle of the plot is from finer sampling near the 5:1 resonance. Center, Right: Keplerian semi-major axis and eccentricity (measured from 
the Pluto-Charon barycenter). The oscillation of these elements demonstrates the non-Keplerian nature of orbits about the Pluto-Charon 
binary. 

B. CODE DETAILS 

We use the Radau integrator (RA15) included with Swifter because it provides a good combination of simphcity, 
accuracy and efficiency. Standard symplectic integrators are not appropriate for this problem because the y assume a 
dominant central mass . Modified symplectic integrators that handle an inner binary have been developed (Chambers 
et al.|2QQ2[ rBeust|2QQ3| ). However for this study we use non-symplectic mtegrators that make no assumptions about the 
masses and orbital architecture of the system. We also tested the standard Bulirsch-Stoer (BS) integrator, another non- 
symplectic method. We opted for the higher order RA15 method because we do not need to handle close encounters. 
Once orbit crossing occurs the system is deemed unstable. 

The main problem we encountered using non-symplectic methods (RA15 and BS) is that they are both inherently 
variable timestep algorithms which the Swifter wrapper occasionally forces to use a fixed timestep in order to produce 
output at regular intervals. (Mercury behaves similarly.) Energy conservation suffers if these forced timesteps are too 
frequent. We were able to mitigate this problem with suitable choices of the output interval, the user-provided dt, 
which we measure in Pluto-Charon (PC) orbital periods. 

If dt is very short compared to the optimal timestep determined by the RA15 algorithm, energy conservation is 
good, but the algorithm's efficiency suffers. At the other extreme, if dt is much longer than optimum RA15 timestep, 
then forced non-optimal timesteps are rare, and energy conservation suffers little. The main cost is that one does 
not have finely sampled output for analysis. For intermediate dt^ many non-optimal timesteps can degrade energy 
conservation significantly. Quantitatively, both long {dt = 1000 PC periods) and short {dt < 0.01 PC periods) 
output times conserved energy to AE/E ^ 10~^ over 10^ PC periods. For comparison, energy conservation drops to 
AE/E ^ 10-^ for dt = 0.1 PC periods. 

For all long integrations, we use a timestep of dt = 1000 PC periods. We set the user-supplied error tolerance 
parameter to 10~^^ for all integrations. The energy error accumulation appears to be, at worst, linear in time. 
Extrapolating to our longest runs, 10^ PC orbits, we expect a total energy error of roughly 1 part in 10^. 

We hope that our description of this issue is of some use for those using non-symplectic codes, especially for long 
timescale integrations where energy conservation is desired. 
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